Coarsening scenarios in unstable crystal growth 
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Crystal surfaces may undergo thermodynamical as well kinetic, out-of-equilibrium instabilities. 
We consider the case of mound and pyramid formation, a common phenomenon in crystal growth 
and a long-standing problem in the field of pattern formation and coarsening dynamics. We are 
finally able to attack the problem analitycally and get rigorous results. Three dynamical scenarios 
are possible: perpetual coarsening, interrupted coarsening, and no coarsening. In the perpetual 
coarsening scenario, mound size increases in time as L ~ f", where the coasening exponent is 
n = 1/3 when faceting occurs, otherwise n = 1/4. 

PACS numbers: 68.55.-a, 05.70.Ln, 81.10.Aj 



Introduction. — The still ongoing large interest towards 
instabilities in crystal growth [l| can be motivated by 
its lying in between the practical problem of controlling 
the growth processes and the fundamental domain of out- 
of-equilibrium transitions and pattern formation. Even if 
some phenomena may be studied in one spatial dimension 
0, most interesting physics of surface growth occurs in 
two spatial dimensions, where rigorous approaches are 
rare and numerical simulations are more problematic and 
much more time consuming. 

The unstable growth of a crystal surface has some simi- 
larities with phase ordering processes, for which the anal- 
ysis of Bray ^ is the most comprehensive and successful 
theory (see also Refs. 0,0] for more recent results). How- 
ever, in two dimensions growth has peculiar features Q, 
because the natural order parameter, the surface slope m, 
obeys a special constraint (V x m = 0), so that result- 
ing domain walls must be straight, which is not the case 
with other phase ordering phenomena. As a consequence 
of that, rigorous analytical approaches to crystal growth 
instabilities are very rare. A paper by Watson and Nor- 
ris 01 , based on the principle of maximal dissipation and 
a couple of papers Q providing exact inequalities for 
the coarsening exponent are some of the few examples. 
Other interesting, albeit heuristic, approaches had con- 
tributed significantly to our understanding of coarsen- 
ing [IMl- 



At the general level, we think that mound formation 
in two dimensional crystal growth, in spite of being quite 
an old problem [ll] , still lacks a systematic and rigorous 
approach. In this manuscript we propose to fill this gap. 
We start by briefly reminding the basic phenomenology 
of unstable crystal growth which is of interest to us and 
presenting a well established continuum description of 
the growing surface. Then the central part of the paper 
follows: the application of the phase diffusion concept 
that allows to face at once a large class of models and 



symmetries. In this approach, the pattern instability is 
signaled by a negative phase diffusion coefficient, which 
implies coarsening, i.e. an increase of the size L of the 
pattern. In the case of perpetual coarsening, the depen- 
dence of the diffusion coefficient on A allows to determine 
the growth law, L(t). Our work gives the first rigorous, 
general framework to study unstable crystal growth and 
we provide important results concerning coarsening sce- 
narios and coarsening exponents. 

Phenomenology and model. — Here we focus on depo- 
sition processes on a high symmetry substrate, where 
atoms or molecules arrive ballistically and thermally dif- 
fuse until they are incorporated into the crystal [l|. A 
key process, which is the possible cause of the instability, 
is the attachment to steps: if adatoms stick to ascending 
steps more efficiently than to descending steps, an up- 
hill current forms and the growing surface destabilizes, 
forming a mound structure. 

Experimental results are available for many systems, 
specially metals, and different dynamical scenarios have 
appeared. Cu(lOO) 14 1 and Rh(lll) [l^ are examples 
of a coarsening process which occurs dur ing the whole 
experimental time scale. Instead, Pt(lll) [16[ is a proto- 
type for coarsening which immediately stops or does not 
even start, while the height of mounds goes on increas- 
ing. In some cases, slopes of mound tend to constant 
values, corresponding to faceting. In other cases, ob- 
served slopes increase as well. Finally, the symmetry of 
the pattern is related to the unstable orientation: (100) 
surfaces produce a 4-fold pattern, while (111) surfaces 



produce a 3- fold pattern 17|. The quantitative determi 



nation of the coarsening exponent, when appropriate, is 
made extremely difficult by the short experimental time 
scale over which L(t) is measured and by the relatively 
smallness of the coarsening exponent, n < ^. 

Theoretical work has accompanied experiments with 
Kinetic Monte Carlo simulations llSi and with continuum 
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approaches [6|, where the dynamics of the local height 
of the surface, /i(x, i), is governed by a partial differ- 
ential equation. Much effort has been devoted to give 
robustness to the equations, while the equations them- 
selves have been mainly studied numerically. Here below 
we propose a more systematic and rigorous approach in 
the continuum limit, starting from the well established 
class of equations: 



dh{x,t) 
dt 



= -V- [j(V/i) + V(V2/i)] = -V-Jt 



(1) 



whose conservative form signals that, at the experimen- 
tal conditions usually set for these deposition processes, 
overhangs and evaporation are negligible. The total cur- 
rent, 3 tot, is composed of two terms: the nonequilibrium, 
slope dependent current, j(V/i), describes the effects of 
the Ehrlich-Schwoebel barrier, which hinders adatoms to 
stick to the lower step, with j ~ V/i for |V/i| <C 1 [l3l |: 
the second term, V(V^/i), is the MuUins term, describing 
termal detachment of adatoms from steps and nonequi- 
librium effects [l^. In Eq. ([T]), the white noise due to the 
flux can be omitted Wc would like to highlight that 
we do not give explicit form to the slope dependent j: our 
analytical approach embraces any current j . We can even 
add that Eq. ^ is also relevant for thermodynamically 
unstable surfaces [20j . 

The method. — Equation ([T]) admits the trivial solution 
h = 0, corresponding to the flat profile. Setting h = 
S exp(a;t-|-ik-x) and linearizing in 5 yields u;{k) = k'^ — k'^, 
k = |k|, showing linear exponential growth, until the 
nonlincarity can no longer be neglected. 

The amplitude growth of mounds is thus fast, while 
wavelength rearrangement (if any) is slow and follows a 
diffusion process, as seen below. At short time scale the 
structure assumes a periodic pattern defined by two ba- 
sis wave vectors qi,q2, or, similarly, by the two phases 
fi = Qi ■ X. The slow evolution of the phase makes it le- 
gitimate to introduce, besides the fast variable x, a slow 
dependence on time and space, T and X, so that the local 
wave vector q = q(T, X). Owing to the diffusion charac- 
ter of the phase, we expect T = e^t and X = ex, where e 
is a small parameter measuring the long wavelength mod- 
ulation (small gradient) of the phase. It is convenient to 
introduce slow phases, related to the fast one by ipi = £(^i, 
so that c{i = Vx(/3i = Vx'0i- ^ multiscale spirit, from 
these definition together with h = ho + ehi + . . . , the 
following expansions for the operators hold: 

dt = e {idTi^i)d^, -f (9tV'2)9^J (2) 
V = Vo+eVx, (3) 

with Vo = qi9^i + and Vx = (9^, 9y)- Finally, 

the current j reads: 



j(V/i) - j(Vo/io) + eJiVahi + Vx/io), 



(4) 



where J7 is the Jacobian matrix. Reporting the above ex- 
pansions into Eq. ([T]) we obtain successively higher order 



contributions in powers of e. To order we gain 

= Vo • [j(Vo/io) + Vo(Vg/io)] = Vo • (Jo)*ot = M[hol 

(5) 

which is a nonlinear equation satisfied by the steady-state 
profile /iQ. Since we are interested in growth on high 
symmetry substrates, the stronger condition (Jo)tot = 
is fulfilled. Then, to first order we obtain: 



C[hi] = g(/io,i/'i,V'2), 

where 

C[hi] = -Vo • [J{Vohi) + Vo(V2^i)] 
is the Frechet derivative of N , and 
g EE {dTiJi)d^,ho + {dTilJ2)d^^ho + Vo • [J(Vx/io) 



(6) 



(7) 



(8) 



+ Vx(V2/io)-HVo(V?/io)]. 
In order to avoid secular terms, we use the Fredholm's 



alternative theorem [2l|, according to which Eq. ([6]) has 
solutions if and only if {v,g) — ]2M, where C^v] = 0. 

While the derivation of the phase equation can be 
achieved for any nonlinear equation, we focus here on 
Eqs. (H]) such that £ is a self-adjoint operator, = C. 
This condition is equivalent to saying that the Jacobian 
i7 is symmetric, which, in fact, is not a restriction: all 
surface currents j discussed in the literature provide a 
symmetric JT". Thanks to the translational invariance of 
J\f with respect to the space variable, we easily obtain non 
trivial solutions of the form Vi = d^pMo, i — 1,2 (Gold- 
stone mode). This leads us to two diffusion equations 
(i = l,2): 



a,/3,7 = l,2 



(9) 



"-7 



where the diffusion coefficients have the following expres- 
sions: 



f>lQ 



(/ll,/ll)(/l2,/l2) - {hi,h2) 



(10) 



and D'p" ^t^^ ^jsj- Here above we have introduced the 
more compact notation hj = d^p-ho = djh^, and 

dho 



dho 



2qjyqil3dadidjho 
S/jyVldaho, 



(11) 



where the use of the zeroth order relation j(Voft.o) = 
— Vo(VoAo) has allowed to replace the current with 
derivatives of the steady profile. 

We conclude this part stressing that our e— expansion 
might continue at higher order, leading to nonlinear 



terms in Eqs. ([9]). We refer the reader to the discussion 
in 0], Sect. VIC. 

Pattern symmetries and phase dijfusion equations — 
For definiteness we focus in the following on square (4- 
fold) and hexagonal (6-fold) patterns. In the square 
case, we set qi = q{l,0) and q2 = 9(0,1) and exploit- 
ing parity and 4- fold symmetry for /ig, physical prop- 
erties are equivalent upon the changes ipi o ip2- In 
the hexagonal case, we have qi = q{l/2,\/3/2) and 
q2 = g(l/2, — -\/3/2) and parity and 6-fold symmetry 
imply invariance under the changes {ipi O ip2) and 
{(fi — > — v?2, V2 </'i + ^2)- The interesting result is 
that, for both symmetries, Eqs. ^ assume a much sim- 
pler form |23| : 



drlpl = ('0l)ll£'ll + (V'l)22£'22 + ii'2)l2Di2 (12a) 
dTi'2 = (V'l)l2£'l2 + (V^2)ll£'22 + (^2)22£'ll, (12b) 

with only three non vanishing diffusion coefficients. The 
hexagonal pattern enjoys an additional property: 

Dn-D22-Di2 = 0, (13) 
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FIG. 1. The amplitude of stationary solutions as function 
of for square symmetry (2gMAx = ^max = 1): the full 
line indicates the branch along which the system will coarsen; 
along the dashed line there is no coarsening. Branches marked 
by the thinnest line correspond to a larger value of C4. Inset: 
the three possible scenarios for the growth model in 2D found 
in the limit k — >■ femAx, in which the sign of the diffusion 
coefficients that can give instability is linked to the behaviour 
of the ampUtude a with respect to the wavelength A. We 
have coarsening (full hne), no coarsening (dashed line), and 
interrupted coarsening (dotted line). 



leaving us with only two independent diffusion coeffi- 
cients. 

Phase instability means wavelength rearrangement, 
leading to coarsening. To investigate stability of the 
pattern we set V'i,2(X, T) = V'i°2 6xp (riT) exp (zK • X), 
where instability is signaled by at least one positive eigen- 
value rii_2(K). In the hexagonal case we straightfor- 
wardly find from Eqs. ([T^ that: 



Q7 



(14) 



61 = 7r/4, for which n^^^iK) = 
and n2^\K) = -{Dn + D22 



-{Dn + D22-Di2)KV2 
- Di2)K^/2. Again, the 



stability is dictated by the sign of diffusion coefficients: 

are negative, while 
have no fixed sign: 



we find that the eigenvalues Q^'^^^ 



the eigenvalues il^'^^^ 



Dii = 



1 



{hi) 



[d,{qHh\,))+q^d,{h\^)+q^{hl^)], (16a) 



Since Dij are expressed in terms of integrals involving the 
steady periodic pattern (see Eq. (fTO|) ). stability can be 
linked to the property of steady-state solutions only. It is 
a simple matter to show that D22 is positive, thus VLi < 0. 
The analysis of Dn is much more involved. After several 
manipulations we find: 



Dll = ^d,{q^/\hl^)) 



ficd^il). (15) 



The sign of Dn (and thus of VI2) is given by the sign of 
the slope of the function A = (9^''l(^i2)) with respect 
to q. In the ID models studied in [4I this function was 
found to be the amplitude of the steady state solutions. 
The present study reveals a new class of dynamics where 
the relevant quantity is wA(g), having a more abstract 
meaning than just the amplitude itself. At present a 
simple physical interpretation of this function is missing. 

For square symmetry, the spectrum is anisotropic, de- 
pending on the angle 9 between K and the X axis. 
However, symmetry imposes that the growth rate may 
be maximal along two directions: (i) = 0, for which 
n\{K) = -D22K^ and n°{K) = -DuK^, and (ii) 



Dii+D 



22 



Dr2 = 



(hi) 



2q'{hl,) 



(16b) 



Coarsening scenarios and coarsening exponents. — 
Even if an exact determination of diffusion coefficients 
can only be done numerically, there are two crucial limits 
allowing analytical treatment: the small amplitude limit 
{k — )■ few AX = 1) and the large wavelenght /large ampli- 
tude limit. The former provides the proof of the existence 
of different coarsening scenarios, the latter provides the 
explicit values of coarsening exponents. 

For the small amplitude limit, we have considered an 
isotropic current of the forrnUm, C2, C4) = m(l -|-C2m-^ + 
C4m*) and we have found 2J] steady solutions and dif- 
fusion coefficients for both square and hexagonal sym- 
metries. In both cases, the phase stability is controlled 
by the amplitude ai{q) of the steady state: we have in- 
stability, i.e. coarsening, if a[{q) < (see Fig. 1). The 
sign of this derivative is determined by the coefficients C2 
and C4: perpetual coarsening occurs if C2 < 0,C4 < (or 
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C2 > 0,C4 < 0), no coarsening if C2 > 0, C4 > 0, coarsen- 
ing up to a critical wavelenght (interrupted coarsening) 
if C2 < 0, C4 > 0. The basic three scenarios are illustrated 
in the inset of Fig. 1. 

Perpetual coarsening is determined by a phase insta- 
bility occurring for any q. In this case we can extract the 
coarsening exponent n, L{t) « t", through the following 
dimensional law: 



1^(9)1 



L2 
t 



(17) 



with q ~ 2Tr/L. The application of Eq. ([T7)) has always 
given values in agreement with known results, both in 
ID and in 2D |g. 

We insist on highlighting that the form of j is arbi- 
trary and that here we shall only distinguish between two 
broad classes of systems: those which exhibit slope selec- 
tion and those which do not. In the first case the current 
j(m) has zeros for finite values of the slope m* (as, for ex- 
ample, with j = m(l — m^)), so that pyramids will grow 
in size but with a slope tending to the constant value m* , 
leading to faceting. This implies that m* = {dxhQjdyho) 
is constant, except along domain walls, that have a fi- 
nite but small thickness. We first consider the square 
symmetry, so, for example: 



Conclusions and discussion — Time and length scales 
of the emerging structure have been investigated within 
the phase diffusion notion, which has allowed to relate 
the growth dynamics to the properties of stationary so- 
lutions. The sign of the coefficients that dictate insta- 
bility is related to some abstract function (see Eq. 
that depends on steady-state solutions. In the limit of 
small amplitude the abstract function coincides with the 
amplitude of the pattern. In this case three kinds of sce- 
narios, namely perpetual coarsening, no coarsening, and 
interrupted coarsening, are possible. This is the first ana- 
lytical evidence of such dynamical scenarios in 2D growth 
models. 

As for the coarsening exponents, they deserve a special 
discussion. In fact, while the result for a faceted 6- 
fold pattern and the result (fT9)) for an unfaceted pattern 



11, 2 



= ciiX^+oiX^) 



where cn is a positive constant. With the same reason- 
ing, at the leading order in A we find = ci2A'^ and 



agree with numerics and scaling considerations 
the exponent n = 1/3 for a faceted square pattern seems 
to be in contradiction with some numerics (25I [26| , where 
a slower coarsening n = 1/4 is reported. These authors 
argue that in the square case the coarsenin g dy namics 
may be slaved to the appearance of roof-tops [26|, a type 
of domain wall which is not present in a regular lattice 
of square pyramids, where only the so called pyramid 
edges appear This type of dislocations, according 

to Hi Q, pla" a major role in their numerical simu- 
lations, leading to a slowing of coarsening: n = 1/4 in- 
stead of n = 1/3. However, more recent results [l^ show 
that such topological defects may or may not appear, 
depending on the explicit form of the slope dependent 
current. Therefore, from the point of view of numer- 



(^1,2) = c^^^- Therefore, we get Dn = -2gci2/c^, that ical simulations, the faceted square case does not give 
is negative, signaling instability. For the other coefficient 
we have (Dn -I- D22 + D12) = 2g(ci2 — cii)/c^, with 
Cii > C12 for phase instability. Plugging the results for 
the Dij 's into ((T7| we finally obtain: 



1/3 



(18) 



A similar analysis for hexagons allows us to determine 
(fT6|l that has exactly the same behaviour with A as for 
the square symmetry, meaning that coarsening exponent 
n = 1/3 also holds for 6- fold symmetry. 

We have also analyzed the coarsening exponents for 
models without slope selection in the steady pattern: 
when L increases, the slope increases as well. Assuming 
that the mound profile changes only along one direction 
while remains constant along the perpendicular one, we 
can find the asymptotic profile. Thus, for a current that 
behaves asymptotically (large slope) as j(m) ~ l/|m|^, 
/3 > 1, we have l/|m|^ = —d'^m/dx^ . By using its solu- 
tion, we have finally obtained the result: 



t 



1/4 



(19) 



for both symmetries and for any values of /3, as in the ID 
case [3]. 



a unique picture. We should also add that in [25|, |26| 
the square pattern was pointed out as metastable and 
therefore defect generation was required to overcome the 
metastability barrier, but we have shown exactly that the 
square pattern is linearly unstable with respect to phase 
fluctuations. 

A theory of coarsening in the presence of topological 
defects is beyond the reach of the present study. How- 
ever, our results support the idea that coarsening can 
take place via pyramid coalescence, without the assis- 
tance of defect generation. Therefore, since defect occur- 
rence was found to slow down coarsening, our analysis 
might suggest that there always exists a faster channel 
for dynamics, i.e. without the intervention of topological 
defects. It would be an important future line of research 
to see whether the presence of defects or the lack thereof, 
is a fundamental ingredient for some specific equations, 
or rather does it depend on initial conditions, boundary 
conditions, system size, etc. . . 

All these results have been established without evok- 
ing any specific form of the slope dependent current j(m), 
therefore pointing to the existence of a universality class. 
Furthermore, the results are independent of the consid- 
ered symmetries: hexagons and squares. It would be 
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an interesting task for future research to see whether or 
not this conclusion survives for the three other Bravais 
lattices (rectangular, rectangular centered and oblique). 
Finally, we point out that unstable growth on a (111) sur- 
face produces a 3-fold pattern, which does not fall within 
Bravais lattices. It is an open question to understand 
how such pattern can be analyzed within our approach. 

A last word for experimental rcults is appropriate. 
While the large variability of experimental values for 
n [l^ prevents from a quantitative comparison, our 
different dynamical scenarios for coarsening are well ob- 
served [343. 

Acknowledgements — CM thanks CNR for a Short- 
Term Mobility award from their International Exchange 
Program and CNES for financial support. Useful discus- 
sions with Matteo Nicoli are also acknowledged. 



* |sofia.biagi@ ujf-grenoble.fr] 

^ chaouqi.misbah@ujf-gre noble.£r| 
j3aolo.politi@isc.cnr.it 
[1] C. Misbah, O. Pierre-Louis and Y. Saito, Rev. Mod. 

Phys. 82, 981 (2010). 
[2] P. Politi and C. Misbah, Phys. Rev. E 73, 036133 (2006). 
[3] A. J. Bray, Adv. Phys. 43, 357 (1994). 
[4] J. J. Arenzon, A. J. Bray, L. F. Cughandolo and A. Si- 

ciHa, Phys. Rev. Lett. 98, 145701 (2007). 
[5] C. Misbah and P. Politi, Phys. Rev. E 80, R030106 
(2009). 

[6] P. Politi, G. Grenet, A. Marty, A. Ponchet and J. Villain, 

Phys. Rep. 324, 271 (2000). 
[7] S. J. Watson and S. A. Norris, Phys. Rev. Lett. 96, 

176103 (2006). 

[8] R. Kohn and X. Yan, Commun. Pure Appl. Math. 56, 
1549 (2003). 

[9] Bo Li, Nonlinearity 19, 2581 (2006). 
[10] M. Rost, J. Krug, Phys. Rev. E 55, 3952 (1997). 
[11] L. Golubovic, Phys. Rev. Lett. 78, 90 (1997). 



[12] L.-H. Tang, P. Smilauer and D. D. Vvedensky Eur. Phys. 

J. B 2, 409 (1998). 
[13] J. Villain, J. Phys. France 1, 91 (1991). 
[14] J.-K. Zuo and J. F. Wendelken, Phys. Rev. Lett. 78, 2791 

(1997). 

[15] F. Tsui, J. Wellman, C. Uher and R. Clarke, Phys. Rev. 
Lett. 76, 3164 (1996). 

[16] M. Kalff, P. Smilauer, G. Comsa and Th. Michely, Surf. 
Sci. Lett. 426, 447 (1999). 

[17] T. Michely and J. Krug, Islands, mounds, and atoms: 
patterns and processes in crystal growth far from equilib- 
rium (Springer, 2004). 

[18] M. Siegert, M. Plischke, Phys. Rev. E 53, 307 (1996). 

[19] P. Politi and J. Villain, in Surface Diffusion: atomistic 
and collective processes, Ed. M. C. Tringides (Plenum 
Press, 1997). 

[20] J. Stewart and N. Goldenfeld, Phys. Rev. A 46, 6505 

(1992) . F. Liu and H. Metiu, Phys. Rev. B 48, 5808 

(1993) . 

[21] D. Zwillinger, Handbook of differential equations. (Aca- 
demic Press, 1998). 
[22] The scalar product is defined as (/, g) := 
1 



[23; 



(2^)^70 Jo 

Since the equality of mixed partials, it is preferable to 
define new diffusion coefiicients 



D 



Pi 



Pita 

Df^ + D. 



7/3 



/? = 7 



(20) 



and to assume /3 < 7 in the diffusion equations ([9]) 
[24] See Supplemental Material 

for 



|http : // link . aps . org/supplement al/ 



the 



at 

details 



of the small amplitude limit. 
[25] D. Moldovan and L. Golubovic, Phys. Rev. E 61, 6190 
(2000). 

[26] M. Siegert, Phys. Rev. Lett. 81, 5481 (1998). 

[27] A pyramid edge is a domain wall where only one com- 
ponent of the slope changes. A roof top is a domain wall 
where both components of the slope change. 

[28] A. Levandovsky and L. Golubovic, Phys. Rev. B 69, 
241402 (2004); L. Golubovic, A. Levandovsky and D. 
Moldovan, East Asian J. Appl. Math. 1, 297 (2011). 



